Dirac and Weyl Equations on a Lattice as Quantum Cellular Automata 



On 



O 



Iwo Bialynicki-Birula 

Centrum Fizyki Teoretycznej PAN* 
Lotnikow 32/46, 02-668 Warsaw, Poland 
and 

Institut fur Theoretische Physik, Johann Wolfgang Goethe-Umversitat 
Robert-Mayer-Strasse 8-10, Frankfurt am, Main, Germany 

A discretized time evolution of the wave function for a Dirac particle on a cubic lattice is 
represented by a very simple quantum cellular automaton. In each evolution step the updated value 
of the wave function at a given site depends only on the values at the nearest sites, the evolution is 
unitary and preserves chiral symmetry. Moreover, it is shown that the relationship between Dirac 
particles and cellular automata operating on two component objects on a lattice is indeed very close. 

£-fj , Every local and unitary automaton on a cubic lattice, under some natural assumptions, leads in the 

^\ • continuum limit to the Weyl equation. The sum over histories is evaluated and its connection with 

^\ ' path integrals and theories of fermions on a lattice is outlined. 
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Oh! 

The aim of this paper is to present a very simple lattice algorithm that reproduces in the continuum limit the prop- 
agation of massless or massive Dirac particles. The discretized time evolution described here satisfies the fundamental 
physical requirements of locality (the updated value of the wave function at a given site depends only on the values 
at the neighboring sites), unitarity (the norm of the wave function is preserved), and chiral symmetry (independent 
_h \ propagation of two hclicity states). 

>• The main, and totally unexpected, result of this investigation is a discovery that for two-component wave functions 

on a cubic lattice the Weyl equation necessarily follows in the continuum limit from locality, unitarity, and two 
additional, natural assumptions: (A) wave functions that have constant values throughout the lattice do not change 
in time and (B) the evolution algorithm preserves the symmetry of the lattice. Thus, the rotation group and spin 
■ emerge automatically in the continuum limit from unitary dynamics on a cubic lattice. 
(*C) I The results presented in this paper are directly related to numerous proposals of path integrals for a Dirac particle 
ON ■ since an iteration of a discretized time evolution gives a sum over histories. The path integral for the Dirac particle 
j in one space dimension was found by Feynman jl],^] and independently by RiazanovJ0|. Even this relatively simple 
problem, where there is no spin to complicate matters, is still attracting attention There is no consensus at all 

O i| as to the form of path integrals for a Dirac particle in three dimensions. Previously proposed path integrals fall into 
three broad categories. The first category PJ^-p^[| comprises those approaches that work by "reverse engineering" 
introducing from the outset the Dirac matrices to describe the spin degrees of freedom. These approaches are often 
combined with the proper-time representation of the Dirac operator. Into the second category |l4|-|ll| fall those 
formulations that derive spin from continuously parametrized space of states related to the rotation group. To the 
third category fll7|-p0| belong all approaches based on anticommuting Grassmannian variables. The sum over histories 
that is obtained from my lattice algorithm is distinct from all these path integrals since the evolution is fully discretized 
and the rotation group emerges only in the continuum limit. 

There is an obvious connection between this work and theories of fermions on a lattice but there is also an essential 
difference in the choice of objects to be calculated. I make no attempt to write down discretized versions of the 
Hamiltonian or the Lagrangian working instead with a discretized version of the evolution operator. In my approach 
common difficulties (breaking of chiral symmetry, doubling of fcrmion species, nonlocality) encountered in formulating 
the dynamics of a Dirac particle on a lattice pl|-p7t and the no-go theorem concerning Weyl particles on a lattice |2f| 
are avoided. 

Algorithms for cellular automata describing quantum system have been studied by Groessing and Zeilinger p9| , 
but their discretized time evolution does not conserve probability. In a recent paper Kostin p0| ] has introduced a 
cellular automaton for the Dirac equation that does conserve probability but his algorithm is nonlinear and it gives 
a linear evolution equation only in the continuum limit. Following the example set by these authors I use the term 
quantum cellular automaton to denote discretized evolution of complex wave functions. I believe that the name cellular 
automaton is justified despite the fact that one of the eight properties of cellular automata required by Wolfram 
does not hold: the states are described by continuous, not by discrete variables. However, my algorithm for updating 
the state of the system is synchronous, homogeneous, discrete in space and in time, deterministic, and spatially and 
temporally local. In my case the name quantum is fully justified by the unitarity of time evolution. 
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Let me begin with a massless particle on a cubic lattice. In my quantum cellular automaton the two-component 
wave function cj)(i, j, k,t) is updated for each time increment At according to the following local algorithm 

tf>{i,j, k,t + At) = W +++ 4>{i + l,j + l,* + l s t) 
+W r ++ _0(i + l,i + l,A-l,t) + --- 

+W___0(i-l,i-l,*;-l ) t), (1) 

where the integers i, j and fc are the coordinates of lattice sites in units of the cell size a (r = (ai, aj, ak)) and all eight 
Ws are 2x2 matrices. This algorithm can be written in a compact form 

0(r,t + At) = £V(h)#r + h,t), (2) 

h 

where the summation extends over a set of vectors h that point from a given site to eight nearest sites. Their 
components in units of a are ±1, 

{h} = {{1,1,1}, {1,1,-1}, {1,-1,1}, {1,-1,-1} 

{-1,1,1}, {-1,1,-1}, {-1,-1,1}, {-1,-1,-1}}. (3) 

Upon evaluating the norm of the updated wave function, one finds that the matrices 17(h) must satisfy the following 
algebraic relations to guarantee the unitarity of the transformation (|^) 

£Vt(h)^(h) = 1, (4) 

h 

(h)W (h + h' - h") = 0, (5) 

h 

where h' and h" ^ h' are arbitrary vectors belonging to the set and the sum in (||) extends only over vectors h 
such that the vectors h + h' — h" are also members of the set. It follows from these unitarity conditions that the 
inverse of (|^) has a similar local form so that my automaton is fully reversible, 

0(r,i) =^W t (h)^(r + h,i + Af). (6) 

h 

Since the inverse of a unitary transformation is also unitary, the hermitian conjugate matrices W' must obey the 
same conditions as the matrices W, 

XV(h)Wt(h)=l, (7) 

h 

53w(h)W^(h + h'-h") = 0. (8) 

h 

Depending on the mutual orientation of the vectors h' and h", the conditions (pi) and (^) have one, two, or four 
terms, as exemplified below 

Wj ++ W___ = 0, (9a) 
Wl ++ W+— + Wl ++ W— = 0, (9b) 



+W]__ + W— = 0. (9c) 



In total there are 8 conditions of the type (|9a|), 12 conditions of the type (pb|), 6 conditions of the type (pq), and one 
condition (Q) and then the same number of conditions with W and interchanged. Nonetheless they can all be 
satisfied by the following simple and symmetric choice of Ws 

W+++ = q + P u W++- = g_P 2 , W+-+ = q_P lt 
W+— = q+P 2 , = q-P 3 , W-+- = q+Pi, 

W—+ = q+P 3 , W— = q-P A , (10) 
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where 



g+ = (l + *)/4, q- = (1-0/4, (11) 

and 

Pi = [ ] °A, p 2 f0 1 
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In order to find the continuum limit of the evolution equation I shall use the exponential representation of the shift 
operators to write (^|) in the form 

<j>(r, t + At) = Y^ W(h) exp(h • V)^(r, t). (13) 

h 

Expanding both sides of this equation in powers of At and a (h is of the order of a) and keeping only linear terms, I 
obtain the Weyl equation, 

d t <f>(r,t)=ca-V4>(v,t), (14) 

where c = a/ At. Incidentally, to obtain a discretized Weyl equation in two space dimensions one may choose the four 

matrices W++, W_| , W |_, and W as equal to the matrices Pi. 

I shall show now that for every set of 2 x 2 matrices satisfying the unitarity conditions, and not just for the particular 
choice (|l^), one obtains the Weyl equation in the continuum limit. To prove this let me first expand the exponential 
operator in (|l3|) into powers of a, 



^iy(h)exp(h- V) = 1 + OS-V+---, (15) 

h 

where I have made use of the assumption A that a homogeneous wave function does not change in time, 

£V(h) = l, (16) 

h 

and I have introduced the matrices s^, 

= w +++ + w++^ + w<v+ + w+— 

(17a) 
(17b) 
(17c) 

Since the differentiations are antihermitian operators, the unitarity of the operator ( |l5[ ) requires that the three 
matrices s, must be hermitian and, therefore, the formulas ( |l7| ) must also hold with all the matrices W replaced by 
their hermitian conjugates. This fact will enable me to use the unitarity conditions (||) and (||) to calculate products of 
the matrices Sj, (i = x, y, z). The unitarity conditions also imply (to see this, one has to write them down explicitly) 
that all products of a IF matrix by its hermitian conjugate commute, and therefore they can be simultaneously 
diagonalized. From the equivalence of all six lattice directions (assumption B) I conclude that the eigenvalues of all 
matrices FF^(h)lF(h) must be the same and then from equations (|J) and (^|) I find that these eigenvalues are 1/4 and 
0. It is now a matter of tedious but otherwise straightforward algebraic manipulations to show that the matrices Sj 
satisfy the familiar anticommutation relations, 
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+ s j s i = 25 ij . (18) 



There exist only two inequivalent two-dimensional representations of these relations: s = <xors = — er. They describe 
the propagation of two helicities. Thus, up to a choice of helicity, the universality of the Weyl equation under the 
listed assumptions is established. 
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For massless particles each helicity state propagates independently. Chiral invariance (or CP symmetry) is expressed 
in my discretized form of time evolution by the fact that the matrices W corresponding to the two helicities are related 
by the spatial reflection 

W(h) -> W(-h), (19) 

as seen from the formulas (|l7|). For massive particles, however, the two helicity states are mixed by the mass term. 
Therefore, the discretized time evolution for a massive particle must be described, as in the standard Dirac equation, 
in terms of two two-component wave functions. The discrete time-evolution algorithm for a massive Dirac particle 
can again be written in the same general form (jfy as for the massless case 

^(r,t + At)=5^D(h)V(r + h,t). (20) 

h 

The 4x4 matrices D that act on four-component wave functions ^(r, t) can be expressed in terms of the matrices 
W(h) and W(— h) by the following block-matrix formulas (from now on I set c = At/a — 1 and % = 1) 

= / cos(ma) idn(ma) W(h) \ (21) 
v wsin(ma) cos(ma) J\ W{— h) J 

One may check that such a modification does not affect the unitarity conditions; the matrices D will also satisfy all 
of them. The continuum limit of (|2^) gives the Dirac equation in the Weyl representation of the Dirac matrices 

d t ip(r, t) = (p 3 a- V + ipim)il>(r, t). (22) 

I shall now write down the sum over histories for the Weyl particle arising from my time-evolution algorithm (|^) 
and I shall directly show that it yields the propagator in the continuum limit. The iV-fold iteration of the single-step 
evolution leads to the formula 

4>(r, t + NAt) 

= W ( h i) • • -W(h N )(f>(r + hi + ■ • • + h N ,t). (23) 

hi h N 

When the sum in this formula is restricted to only those combinations of the vectors that produce a given displace- 
ment, we will obtain a discrete version of the propagator 

K(t - r', t) = Wthi)- ■ ■W(h N )6 ry+bl+ ... +ilN , (24) 

hi,...,h N 

where t = N At. Each term in this sum corresponds to a lattice trajectory that in N steps connects the initial and 
the final lattice sites. A single step described by a vector h is represented by the matrix W(h). Contributions from 
all trajectories are coherently added. 

The constraint on the sum in ( p4| ) can be handled (cf., for example, J3^]) with the help of the Fourier representation 
of the Kronecker delta, leading to the following integral form of the propagator 

K(r-r',t) - (^) 3 /d 3 fcW^(k)exp[jk-(r-r')], (25) 

27T J 

where 

W(k) = Y W(h) exp(ik-h), (26) 

h 

and the integration extends over the Brillouin zone: — ir/a < ki < ir/a. 

The integral (^5|) can be evaluated with the use of the explicit representation (|l0|) of the matrices W, 

VF(k) — CxCyCz ~i~ S'x^ySz ^\(,SxCyCz Cx$y Sz)@x 

~^~(.CxSyCz ~t~ SxCySz)@y ~t~ (CxCy^z $x Sy Cz )&z\ i (^^) 

where 
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Cj = cos(fcja), Si — sin(feja). 



(28) 



This is a unitary matrix and its eigenvalues X± are 

A± = exp(±z<p) 

— C-xCyCz ~f~ SxSy^z ^ (QrCyC^ ~b Sx^ySz*)^- (^^) 

An easy way to evaluate the iV-th power (JV = t/At) of the matrix (p7|), needed to calculate the propagator, is to 
write W in terms of its eigenvalues and its unitary diagonalizing matrix V, 

W(k)= V ( eM ' ip) (3°) 
\_ exp(— up) J 



exp(— itp/a) 
In the limit, when a — > 0, one obtains 

lim(<p/a) = |k| =: k, (32) 



and finally 



lim V = (ker - ka z )/J2k(k - k z ), (33) 

a— »0 



(W(k))st = exp(ik-o-t), (34) 



in accordance with the Weyl equation. Thus, the sum over histories ( |23| ) reproduces correctly the propagator in the 
continuum limit. 

Having reproduced the propagation of a massless particle I can easily include the mass in the sum over histories 
since that part has been already solved by Feynman |j]j|] in one dimension. More recently Feynman's "checkerboard" 
picture of a particle zigzagging through spacetime, reversing its helicity at each bend, has also been described by a 
Poisson process [p| Jl5| , pO| , |32] | . This Poisson process must be combined with the propagation of definite helicity states. 
Therefore, the propagator for a massive Dirac particle will be a sum of terms, each term describing a fixed number 
of helicity reversals. Between the reversals the propagation is described by the sum over histories (^4|) separately for 
each helicity. 

Interaction with the electromagnetic field may be accounted for in the standard fashion by changing the phase of 
the wave function at each propagation step according to the local value of the electromagnetic potential. 

I would like to thank Jerzy Kijowski and Joachim Rcinhardt for discussions and Walter Greiner for his hospitality 
at the University of Frankfurt. 
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